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Abstract. The Mock LISA Data Challenge is a worldwide effort to solve the 
LISA data analysis problem. We present here our results for the Massive Black 
Hole Binary (BBH) section of Round 1. Our results cover Challenge 1.2.1, where 
the coalescence of the binary is seen, and Challenge 1.2.2, where the coalescence 
occurs after the simulated observational period. The data stream is composed 
of Gaussian instrumental noise plus an unknown BBH waveform. Our search 
algorithm is based on a variant of the Markov Chain Monte Carlo method that uses 
Metropolis-Hastings sampling and thermostated frequency annealing. We present 
results from the training data sets and the blind data sets. We demonstrate that 
our algorithm is able to rapidly locate the sources, accurately recover the source 
parameters, and provide error estimates for the recovered parameters. 
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1. Introduction 

Massive black hole binaries (BBH) are expected to be one of the strongest candidate 
sources for the Laser Interferometer Space Antenna, LISA, a joint ESA-NASA 
mission that will search for gravitational waves (GW) is the frequency bandwidth 
10~ 5 < //Hz < 1 pQ. The detection of BBHs by LISA is important for many reasons. 
Firstly, it will allow us to carry out a test of gravity in the highly nonlinear strong- 
field regime [21 EH]. Secondly, it will allow us, in conjunction with other astronomical 
methods, to investigate such things as galaxy interactions and mergers out to very 
high redshift (z > 10). It will also allow us to test galaxy formation models such as 
Hierarchical formation, where it is believed that modern day galaxies were formed from 
the merger of smaller "seed" galaxies. Due to their high masses, the inspiral phases for 
these systems occur at frequencies which are unavailable to the ground based detectors 
due to their low frequency cutoffs. Also, unlike galactic binaries and Extreme Mass 
Ratio Inspirals (EMRIs), the BBHs are very clean sources with detectable signal to 
noise ratios (SNR) of order ~ 10 — 1000s. This means that we will not have to deal 
with confusion noise between sources [5], something which is very important in the 
search for galactic binaries and to a lesser extent in the search for EMRIs. 

Here we describe our analysis of the BBH component of the Mock LISA Data 
Challenge j6]. Our search algorithm is a variant of the Markov Chain Monte Carlo 
method (MCMC) which has been described in References [3 [HE]. The algorithm 
uses a mixture of frequency annealing with thermostated heat, simulated annealing, 
plus a 5-D exploration of the posterior distributions for the search parameters (we 
use a generalized ^-"-statistic to automatically search over the distance, inclination, 
polarization and initial phase). We refer the reader to Reference [5] for a full 
description of the algorithm. 

In our earlier work we used the Low Frequency Approximation (LFA) [10] to 
model the instrument response, but to our surprise this proved to be insufficiently 
accurate for the MLDC data sets where the full detector response is used. Since the 
maximum frequencies of the injected signals were below 2 mHz, we had expected the 
LFA to be adequate, but when running on the training data sets we found systematic 
offsets in many of the parameters. The parameter recovery improved significantly 
when we upgraded our instrument response to the Rigid Adiabatic Approximation 
(RAA) [11], which includes finite armlength effects. Our interpretation of this finding 
is that while the differences between the LFA and RAA are small at 2 mHz, the 
differences are amplified by the very large contribution to the signal to noise ratio that 
comes from from the final cycles of the inspiral. Indeed, we suspect that the remaining 
small systematic offsets in the recovered masses can be traced to the approximate 
nature of the RAA. 

We modified our barycenter waveforms and signal tapers to agree with those used 
to inject the simulated signals |12j . but rather than searching over initial phase (the 
phase parameter used to generate the waveforms), we continued our earlier practice of 
searching over the phase at coalescence, as we have found this to give better acceptance 
rates in the search chains. 

The organization of the rest of the paper is as follows : In Section [2] we present a 
very short discussion of the search algorithm used. We define the Metropolis-Hastings 
sampling, the frequency annealing with thermostated heat and the simulated annealing 
scheme used. Section [3] contains a presentation of results for both the training and 
blind data sets for the challenges. 
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2. The Search Algorithm. 

Our search algorithm, which is again explained in more detail in [5] is based on a 
Metropolis-Hastings sampling, which is the central engine of the standard Markov 
Chain Monte Carlo method (MCMC). Our method incorporates both simulated and 
frequency annealing schemes at various stages in the search. The algorithm uses a 
number of proposal distributions to jump around the parameter space, as well as a 
maximization over the time to coalescence. To quickly summarize : Starting with the 
signal s(t) and some initial template h(t,x) the (un-normalized) posterior density at 
x is computed. We then draw from a proposal distribution and propose a jump to 
another point in the space y. The posterior and proposal densities at x and y are then 
compared by forming the the Metropolis-Hastings ratio 

H = K(y)p{s\y)q(x\y) 
w(x)p(s\x)q(y\x) 

Here n(x) are the priors of the parameters, p(s\x) is the likelihood, and q(x\y) is the 
proposal distribution. This jump is then accepted with probability a = mm(l, if), 
otherwise the chain stays at the proposal point. 

The likelihood was calculated using a generalized ^"-Statistic [IS] , which allows for 
an analytic extremization over the extrinsic parameters: luminosity distance; orbital 
inclination, polarization and orbital phase at coalescence. Consequently, the Markov 
Chain portion of our search returns a marginalized posterior distribution. For future 
Challenges, we intend to follow the initial detection with a full parameter search so as 
to obtain the full posterior. 

In the first section of the search we use a frequency annealing scheme to speed 
up the search. The search templates and the inner products used to compute the 
JF-Statistic are terminated at a cut-off frequency / cu t, which is initially set at 4 x 10~ 5 
Hz. The cut-off frequency is then increased as the chain progresses according to 

!lQ-B(l-i) f f < f 

iu J max J ^ J max 

(2) 
./max / ^ /max 

where i is the number of steps in the chain, B is a growth parameter and / ma x is the 
maximum frequency the signals could reach give the priors on the masses etc.. As the 
cut-off frequency is incremented, more of the BBH signal is revealed, and the SNR of 
the best fit templates increases. Thus, frequency annealing acts in a similar way to 
traditional simulated annealing, but with the added benefit of saving in the cost of the 
template generation and J-"-Statistic computation. In testing we found that the search 
chains would sometimes lock onto secondary maxima during the frequency annealing 
phase, so we introduced a "thermostating" procedure to control the effective SNR of 
the recovered signals. This was done by multiplying the noise spectral density in the 
noise weighted inner products by a "heat" factor p, which was adjusted based on the 
SNR of the current template: 

1.0 < SNR < 20 

> { , (3) 

^) 2 SNR > 20 
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Once the frequency annealing stage was completed, the chain was cooled using the 
simulated annealing scheme 



with £ = log 10 /3 max where /3 max is the heat factor at the end of the frequency 
annealing stage. The index j counts from the end of the frequency annealing, and 
the cool down lasts N c steps. A standard MCMC exploration of the marginalized 
posterior distribution function (PDF) commences once (3 = 1 and continues for N e 
steps. Finally, over the course of Nf steps, we freeze the chain to a heat of f3 = 0.01 
to aid in the extraction of Maximum Likelihood Estimates (MLEs) for each of the 
parameters. 

3. Conducting the Challenge. 

The MLDC for binary black holes was broken into two parts. Challenge 1.2.1 had a 
source that coalesced inside the observation period, while Challenge 1.2.2 had a source 
coalescing outside of the observation period. The data sets provided consisted of 2 21 
data points sampled every 15 seconds giving approximately one year of data. 

3.1. Challenge 1.2.1. 

Priors on the source parameters were given for each challenge. For Challenge 1.2.1 
we were told that mass ranges were restricted such that mj <E [1,5] x 10 6 M Q and 
m 2 = mi/x, where 1 < x < 4, and that the time to coalescence lay in the range 
t c e [5, 7] months. We were also told that the system would have a signal to noise in 
the range 450 < SNR < 500 in one interferometer. No priors were given on the other 
source parameters. In addition to the blind data set, a training data set was made 
available with parameters drawn from the same set of priors. 

For the 1.2.1 training data set we ran several 50,000 point chains. These chains 
were composed of a 10,000 point frequency annealing search, a 10,000 point simulated 
annealing phase cooldown, a 20,000 point MCMC chain to explore the PDFs, and 
a final 10,000 point freezing of the chain to extract the MLEs. In Figure [T] we plot 
three different search chains for eight of the nine waveform parameters in the 1.2.1 
training set. We have omitted the search chain for the initial phase. We plot the 
extrinsic parameter chains even though their values were determined by analytical 
extremization at each point in the chain. While the angular variables employ a 
linear scale for the number of steps in the search, the other parameters are plotted 
against a logarithmic scaling to highlight the early convergence that occurs for the 
parameters such as the time to coalescence and the masses. The training data set 
provided for Challenge 1.2.1 proved to be more challenging than similar test cases we 
generated for ourselves. In particular, the sky position of the source was not recovered 
to the accuracy predicted by a Fisher matrix estimate of the measurement errors. We 
attribute this to an unfortunate alignment of the source with respect to the LISA at 
the time of coalescence. At coalescence the motion of the detector turns out to be 
perpendicular to line of sight to the source, so there is little or no Doppler shift of the 
gravitational waves right at the time when most of the SNR is accumulating. Thus, 
directional information gets less weighting than is typical, and the small differences 




(4) 
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Figure 1. A plot of three search chains for eight of the nine parameters in the 
Challenge 1.2.1 training set. The solid lines in each cell denote the true value. 



between the RAA and the full response used to generate the data sets is amplified. 
The problem went away when we tested our search algorithm on the same source 
using data that was simulated with the RAA. In future Challenges we plan to use a 
full detector response model to generate the search templates. 

The sky location for the 1.2.1 blind data set proved to be easier to pin down and 
we were able to get away with shorter search chains of 25,000 points. These chains 
were composed of a 5,000 point frequency annealing search, a 5,000 point simulated 
annealing phase cooldown, a 10,000 point MCMC chain to explore the PDFs at unit 
heat, ending with a 5,000 point freezing of the chain to extract the MLEs. We have 
plotted the search chains for the blind data set in Figure [21 along with the values of the 
injected source parameters which were revealed after we had submitted our results. 
We see that the chains typically find the three most important parameters (M c ,fj,,t c ) 
in around 1000 steps. Once again the sky positions took longer to converge, but in 
this instance the chains converged on the injected parameter values. The search also 
accurately recovered the extrinsic parameters, save for the initial phase. The failure to 
recover the initial phase was due to a bug in our ^-Statistic routine that we overlooked 
when working on the training data. We should also mention that while the recovered 
value for the polarization ip is off by ir from what was injected, the polarization angle 
is only defined up to multiples of it so the solution is physically identical. 

In Figure [3] we have plotted the marginalized PDFs for the extrinsic parameters 
based on the merged MCMC chains from each different run. The solid line in each 
cell is the Fisher matrix prediction for the marginalized posteriors. While the chains 
we not long enough to fully characterize the posteriors, we see that the merged chains 
show good agreement with Fisher predictions for (M c ,8,cf>). The Fisher predictions 
for (//, t c ) do not agree very well with the MCMC results. At this time we do not have 
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Figure 2. A plot of three search chains for eight of the nine parameters for the 
blind Challenge 1.2.1 data set. The solid lines in each cell denote the true value. 



an explanation for the disagreement. 

In Table[T]we compare that values of the key file against the MLEs for the training 
set (top) and the blind set (bottom) . We also quote the 1-a error estimation from the 
Fisher matrix, and the error in multiples of the 1-a error estimates from the Fisher 
matrix. In the table, t c is recorded in seconds. While again we did not search for them, 
the MLEs obtained for (i,1nD L ,ip) give errors of (3.86 x 1(T 3 ,6.62 x 1(T 2 , 1.43) and 
(1.7 x 10~ 2 ,4.15 x 10~ 2 , —2.7 x 10~ 4 ) for the training and blind data setsrespectively. 
The ip value is obtained from converting ip — > ip + 7r. For the blind data set, the source 
had a combined SNR of 664.78, while we recovered a SNR of 658.38 . Again the 
mismatch is due, we feel, to a phasing issue in the code which has since been rectified. 
Each run took approximately 24 hours on a single Mac G5 processor. 

3.2. Challenge 1.2.2. 

The priors for Challenge 1.2.2 gave a time to coalescence of 400 ± 40 days and masses 
chosen such that m\ £ [1,5] x 10 6 M Q and again m 2 = mi/x, where 1 < x < 4. No 
priors were given on the other six parameters, save that the source would have an SNR 
in the range of 20 < SNR < 100 in one interferometer. For both the training and 
blind sets, we used 25,000 point search chains. These chains were composed of a 5,000 
point frequency annealing search, a 5,000 point simulated annealing phase cooldown, 
a 10,000 point MCMC chain to explore the PDFs at unit heat, ending with a 5,000 
point freezing of the chain to extract the MLEs. 

In Figure [4] we plot three search chains for the training data set. The chirp mass 
and time to coalescence are recovered in under 1000 steps. For this particular source, 
the chains also rapidly converge on the correct sky position. Again, while not explicitly 
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Figure 3. A comparison of the marginalized PDFs from the MCMC chains for 
the intrinsic parameters of 1.2.1 against the Fisher prediction (solid line) from the 
mean of the chain. The means of the chain have been subtracted and the values 
scaled by the square root of the variances of the chains. 



searched for, the search correctly recovered the inclination of the orbit and luminosity 
distance. However, the recovered polarization angle is not very good. We should note 
here that we have mapped the key file value back into a to 7r range. 

In Figure [5] we plot three search chains for the blind data set. All five intrinsic 
parameters lock-in after a few thousand steps. Except for the luminosity distance, 
the extrinsic parameters are far from the true values. We attribute this to the fact 
the phase at coalescence is essentially undetermined for systems where we do not see 
coalescence. In Figure [5] we plot the marginalized PDFs for the extrinsic parameters 
based on the merged MCMC chains from multiple runs. The solid line in each cell is 
the prediction of the Fisher matrix at the mean of the chain. There is a discrepancy 
between the Fisher prediction and the posteriors from the chains. A visual inspection 
of the chains indicates a large autocorrelation that extends over thousands of points. 
This slow mixing of the chains implies that we would need to run much longer MCMC 
segments in order to get meaningful posterior distributions to compare to the Fisher 
predictions. 

In Tabled] we compare that values of the injected parameters against the MLEs. 
We again quote the l-cr error estimation from the Fisher matrix and the error in 
multiples of the Fisher matrix based on the MLEs. Also, based on the ^"-Statistic 
search for the intrinsic parameters, the MLEs obtained for (o, In Dl,i/j) give errors of 
(-0.0313,0.0124,1.09) and (-0.71,0.273,-2.16) for the training and blind data sets 
respectively. The source had a combined SNR of 106.54, while we recovered a SNR of 
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Table 1. This table compares the injected parameter values, the MLE for each 
parameter, the la error predicted by the Fisher matrix at the injected values, and 
the difference between the MLEs and the injected values in multiples of the Fisher 
la error estimate the Challenge 1.2.1 training set (top) and blind set (bottom). 




Figure 4. A plot of three search chains for eight of the nine parameters of the 
Challenge 1.2.2 training data set. The solid lines in each cell denote the true 
value. 



105.72 . Each run took approximately 6 hours on a single 2 GHz Mac G5 processor. 
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Figure 5. A plot of three search chains for eight of the nine blind parameters in 
Challenge 1.2.2. The solid lines in each cell denote the true value. 
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Table 2. This table shows the injected parameter values, the MLE for each 
parameter, the la error predicted by the Fisher matrix at the injected values 
and the difference between the MLEs and the injected values in multiples of the 
Fisher la error estimate for the Challenge 1.2.2 training set (top) and blind data 
set (bottom). 
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Figure 6. A comparison of the marginalized PDFs from the MCMC chains for 
the intrinsic parameters of 1.2.2 against the Fisher prediction (solid line) from the 
mean of the chain. The means of the chain have been subtracted and the values 
scaled by the square root of the variances of the chains. 



4. Conclusion. 

The Mock LISA Data Analysis Challenge data sets for binary black hole inspirals 
proved to be a valuable testing ground for our Metropolis-Hastings search algorithm. 
Overall the algorithm performed very well, recovering the injected parameters to 
an accuracy largely consistent with the theoretical error margins. There were also 
some surprises, such as needing to go over to a more accurate treatment of the 
instrument response in our template generation, and a problem with our conventions 
when extracting the initial orbital phase. Work is now underway on the much more 
complicated multi-source Challenge 2 data sets, and early results from runs on the 
training data look very promising. 
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